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We investigate small equal-mass two-component Fermi gases under external spherically symmetric 
confinement in which atoms with opposite spins interact through a short-range two-body model po- 
tential. We employ a non-perturbative microscopic framework, the stochastic variational approach, 
and determine the system properties as functions of the interspecies s-wave scattering length a 3 , the 
orbital angular momentum L of the system, and the numbers Ni and N2 of spin-up and spin-down 
atoms (with Ni — N2 = or 1 and N < 6, where N — Ni + N2). At unitarity, we determine the 
energies of the five- and six-particle systems for various ranges ro of the underlying two-body model 
potential and extrapolate to the zero-range limit. These energies serve as benchmark results that 
can be used to validate and assess other numerical approaches. We also present structural properties 
such as the pair distribution function and the radial density. Furthermore, we analyze the one-body 
and two-body density matrices. A measure for the molecular condensate fraction is proposed and 
applied. Our calculations show explicitly that the natural orbitals and the momentum distributions 
of atomic Fermi gases approach those characteristic for a molecular Bose gas if the s-wave scattering 
length a s , a 3 > 0, is sufficiently small. 



PACS numbers: 03.75.Ss,05.30.Fk,34.50.-s 



I. INTRODUCTION 

Over the past few years, the interest in small trapped 
Bose and Fermi gases, and mixtures thereof, has in- 
creased tremendously for a number of reasons. First, 
atomic gases provide an ideal platform for investigating 
phenomena related to Efimov physics [IH3]- While the 
majority of investigations of the Efimov effect have fo- 
cused on the three-body system, larger systems have at- 
tracted considerable attention recently from theoretical 
and experimental groups [44l5j. Second, small trapped 
atomic systems can be realized by loading an atomic gas 
into an optical lattice flB4Tii ] . If the tunneling between 
lattice sites is small and if the interactions between neigh- 
boring sites can be neglected, then each lattice site pro- 
vides a realization of a trapped few-body system. In this 
setting, one interesting prediction is that effective three- 
and higher-body interactions should emerge [2(|. Third, 
small atomic gases can be viewed as a bridge between 
two-body and many-body systems (see, e.g., Refs. [2TI — 
|24|). In most cases, the two-body system is well charac- 
terized, making a bottom-up approach attractive. Such 
an approach treats increasingly larger systems and even- 
tually connects observables for mesoscopic systems with 
those predicted by many-body theories, e.g., through the 
use of the local density approximation. Fourth, few- 
body systems often times allow for highly accurate treat- 
ments, thereby providing much needed benchmark re- 
sults. For example, a number of lattice-based approaches 
are presently being applied to trapped cold atom systems 
(see Refs. [25H30| for lattice-based treatments of the ho- 
mogeneous system). While these approaches promise to 
be very powerful, currently only a few benchmark results 
are available that allow for a careful assessment of their 
validity regimes. 



This paper treats equal-mass two-component Fermi 
gases under external harmonic confinement with short- 
range s-wave interactions. Our work builds on 
the rapidly expanding number of papers that treat 
trapped three-dimensional few-fermion systems (see, e.g., 
Refs. [HQ, EMU)- The ground state of trapped equal- 
mass two-component Fermi gases, e.g., has been inves- 
tigated numeric ally b y the fixed-node diffusion Monte 
Carlo approach | 2l}- -23. 39] and the stochastic variational 
approach [2l|, |22j, |3(j, l4fj, |42| . In the strongly-interacting 
unitary regime, the properties of the system — motivated 
by analytical treatments that exploit the scale invari- 
ance of equal-mass Fermi gases at unitarity [3TI l32l|— 
have been interpreted within the hyperspherical frame- 
work fU [2^]. In some cases, the excitation s pec trum 
at unitarity has also been investigated [2l|, [22I l32l [33| . 
In addition, small two-component Fermi gases have been 
investigated as a function of the s-wave scattering length 
a s [H IMM E3fil. For small \a s \, a s < 0, the en- 
ergy crossover curve has been analyzed by applying first 
order perturbation theory to a weakly-attractive atomic 
Fermi gas [H, [37], 53 ■ For small \a s \, a s > 0, in contrast, 
the energy crossover curve has been analyzed by apply- 
ing first order perturbation theory to a weakly-repulsive 
molecular gas [11 [37|,l!3 (see also Refs. jHfi|). Small 
two-component Fer mi g ases have also provided the first 
high precision tests [40|] of the Tan relations [46M48| that 
apply to both inhomogeneous and homogeneous s-wave 
interacting Fermi gases. 

Following up on our earlier work, this paper presents 
new results for trapped equal-mass Fermi gases with up 
to N = 6, where N = Ni + N 2 and Ni - N 2 = or 
1. Our main results are: (i) We report extrapolated 
zero-range energies for five- and six-particle systems with 
(N±, N2) — (3, 2) and (3, 3) for various angular momenta 
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at unitarity. (ii) We present energy crossover curves for 
the (Ni,N 2 ) — (3,2) system for the ground state and 
various excited states, (iii) We present a detailed anal- 
ysis of the dependence of the few-particle energies on 
the range of the underlying two-body potential, (iv) We 
present structural properties for the (iVijiVb) = (2,1), 
(2,2), (3,2) and (3,3) systems throughout the crossover, 
including unitarity. (v) We quantify the correlations 
of few-fermion systems by analyzing the one- and two- 
body density matrices as well as the momentum distribu- 
tions. In particular, we propose a measure of the molec- 
ular condensate fraction and apply it to few-fermion sys- 
tems with up to N = 6 atoms. Related analyses have 
previously been pursued for bosonic gases p9M51] and 
one-dimensional systems [52- 54], but we are not aware 
of analogous studies for trapped three-dimensional two- 
component Fermi gases. 

Section [TT] introduces the system Hamiltonian and the 
stochastic variational approach employed to solve the 
time-independent Schrodinger equation for small trapped 
two-component systems. In addition, Sec. HT1 reviews the 
definitions of the one- and two-body density matrices 
and their relationship to the natural orbitals and mo- 
mentum distribution. Section [IIII presents and interprets 
our results for various parameter combinations. Lastly, 
Section IIVI summarizes our main results and concludes. 
Mathematical derivations and discussions of technical as- 
pects are collected in Appendices [A] through [C] 



II. THEORETICAL BACKGROUND 



We take the range r$ to be much smaller tha n the har- 
monic oscillator length aho, where a\ m — y/H/(m a u>). 
The depth Vq, Vq > 0, and the range r are adjusted 
so that the free-space two-body s-wave scattering length 
a s takes on the desired value. We restrict ourselves to 
two-body potentials that support no free-space s-wave 
two-body bound state and one free-space s-wave two- 
body bound state for negative a s and positive a s , respec- 
tively. If the scattering length a s is notably larger than 
the range ro, then the properties of small trapped two- 
component Fermi gases are universal, i.e., independent of 
the details of the underlying two-body potential [HI, [55l — 
[Hj]. Thus, we limit ourselves to parameter combinations 
with vq <C a s and r$ -C a\ w . For these parameter com- 
binations, energy shifts due to p-wave or higher partial 
wave scattering between unlike fermions are negligible. 
In a few cases, we perform calculations for different ro 
and explicitly extrapolate to the ro — > limit. 

Our goal is to solve the time-independent Schrodinger 
equation for the Hamiltonian given in Eq. ([T]), and 
to analyze the energy spectrum and structural proper- 
ties. To this end, we use that the total wave function 
V'totCru • • • ,tn) separates into a relative part tp re \ and a 
center-of-mass part tp cm . The relative wave function i/) re i 
is written in terms of Jacobi vectors p\, ■ ■ ■ , pN-xl its de- 
termination through the stochastic variational approach 
is reviewed briefly in the next subsection. Throughout, 
we assume that center-of-mass excitations are absent, i.e., 
we assume that the center of mass wave function is given 
by 



A. System Hamiltonian 



l/'cm(i?cm) = A^cmexp 



Hi 



Ho/*, 



(3) 



Our model Hamiltonian that describes equal-mass two- 
component Fermi gases with N\ spin-up and N2 spin- 
down atoms (N = N1 + N2 and N% > N2) under external 
spherically symmetric harmonic confinement with angu- 
lar trapping frequency oj reads 



3 = 1 



-ti 2 1 

-Vt + -m a w 2 r=f I + 



2m Q r > 2 



Ni N 

E E ^b(r^). (1) 

j=l k=N ± +l 



where iV cm denotes a normalization constant and i? cm the 
center of mass vector, i? cm = Ylj=i ^j/^- The relative 
wave function tp le i is a simultaneous eigen function of 
the relative Hamiltonian H rc \, the square of the relative 
orbital angular momentum operator, the z-projection of 
the relative orbital angular momentum operator and the 
parity operator. Correspondingly, ip Te i and the associated 
eigen energies E re \ are labeled by the quantum numbers 
L, M L and n. 



B. Stochastic variational treatment 



Here, m a denotes the atom mass and fj the position vec- 
tor of the jth particle measured with respect to the trap 
center (with rjk = \rj — fk\)\ the first N% position vectors 
correspond to the spin-up atoms and the last N 2 posi- 
tion vectors to the spin-down atoms. Hamiltonian (TT]) 
assumes that like fermions are non-interacting. The in- 
terspecies interactions are modeled through a purely at- 
tractive Gaussian two-body potential Vth (r) , 



V t b(r) = -Vbexp 



V2: 



ro 



(2) 



To determine the relative eigen functions tl> Te \ and rel- 
ative eigen energies E le \, we employ the stochastic vari- 
ational (SV) approach 65-68|. Our implementation fol- 
lows that described in Refs . |22|, l36l l42j . and here we only 
emphasize a few key points. The SV approach expands 
the relative wave function tp re \ in terms of a basis set. 
The basis functions themselves are not linearly indepen- 
dent, and the determination of the eigen energies requires 
the solution of a generalized eigen value problem that in- 
volves the Hamiltonian matrix and the overlap matrix. 
Just as with other basis set expansion techniques, the 
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SV approach results in a variational upper bound to the 
exact eigen energies, i.e., to the ground state energy and 
to the energies of excited states. For the interaction and 
confining potentials chosen in this work, the functional 
forms of the basis functions allow for an analytical eval- 
uation of the Hamiltonian and overlap matrix elements. 
The proper fermionic symmetry of the basis functions 
is ensured through the application of a permutation op- 
erator A. For the (3,3) system, e.g., A consists of 36 
permutations (6 permutations each are required to anti- 
symmetrize the three spin-up and the three spin-down 
fermions). 

While the functional forms of the basis functions are 
relatively simple, they are sufficiently flexible to describe 
short-range correlations that develop on a length scale 
of the order of the range r$ and long-range correlations 
that develop on a length scale of the order of the oscil- 
lator length aho [Hj]. This is achieved through the use 
of a comparatively large number of variational parame- 
ters that are optimized semi-stochastically for each basis 
function. In this work, we employ basis functions that 
are characterized by N(N-l)/2 to N{ N- l)/2+3(iV- 1) 
parameters [see Eq. (jAlj) of Appendix E] for an explicit 
expression for the basis functions with L n = + sym- 
metry and Eq. (6.27) of Ref. [67|, or Eqs. (36) and (37) 
of Ref. [42| , for an explicit expression of the basis func- 
tions with arbitrary L employed in this work]. Gener- 
ally speaking, the treatment of states with L n = + 
is numerically less challenging than that of states with 
other symmetries. As the range rg decreases or N in- 
creases, the numerical complexity of the calculation in- 
creases. Also, for a given tq and (A^A^) combination, 
the numerical complexity increases with increasing an- 
gular momentum. The largest calculation reported in 
Sec. IHII uses A& = 3000, where Nj, is the number of fully 
anti-symmetrized basis functions. In many cases, how- 
ever, the optimization procedure of the variational pa- 
rameters is more important than the size of the basis set 
itself. In our implementation, e.g., a notable fraction of 
the computational efforts is directed at optimizing the 
basis functions for a given Nf, as opposed to increasing 
Nf,. The motivation for keeping the basis set relatively 
small is two-fold. First, the use of highly optimized basis 
functions mitigates essentially all problems that would 
otherwise arise from the linear dependence of the basis 
functions [36l . Second, the computational time re- 
quired to calculate structural properties increases with 
increasing Nf,. 

To calculate structural properties, we follow two dif- 
ferent approaches. Where possible, we determine the 
matrix elements for a given operator A analytically, 
and determine the quantity (V'tot|^|^tot)/(V , tot|'0tot) by 
simply adding all matrix elements weighted by the ap- 
propriate expansion coefficients. This approach scales 
quadratically with Nb and is, in most cases, more ef- 
ficient than our second approch, a variational Monte 
Carlo calculation [69[ that uses the wave function opti- 
mized by the stochastic variational approach. In partic- 



ular, we calculate structural observables by performing a 
Metropolis walk that samples the probability distribution 
l^tot | 2 / (V'tot | iptot) ■ The expectation value of A is then de- 
termined by averaging over many possible realizations of 
the system. At each step, the density IV'tot] 2 needs to 
be calculated, resulting in a scaling of the computational 
effort with NbN sala p\ e , where the number of Metropo- 
lis steps A^sampio is generally much larger than iV&. Ap- 
pendix [B] details the Monte Carlo sampling scheme for a 
number of observables that quantify the correlations of 
the system. A distinct advantage of the Metropolis sam- 
pling approach is that it allows for the evaluation of "con- 
ditional observables" such as the quantity p re d{R', R), 
defined below Eq. (fT5|) , for which analytical expressions 
of the matrix elements are not available. 



C. Density matrices, occupation numbers and 
momentum distribution 

To quantify the correlations of trapped few-fermion 
systems, we consider the radial and pair distribution 
functions as well as the one- and two-body density matri- 
ces [7014731 ] . The density matrices not only lead to a prac- 
tical route to determine the momentum distributions as- 
sociated with the spin-up and spin-down atoms, but also 
serve to quantify the non-local correlations of the system. 
For example, for trapped single-species Bose gases, an 
eigen value of the one-body density matrix of the order of 
1 signals a large condensate fraction [4jl [7l|, [73|, l74j. The 
situation is different for two-component fermions |72ll73j. 
Because of the anti-symmetric many-body wave function, 
none of the natural orbitals associated with the one- 
body density matrix can be occupied macroscopically. 
In fermionic systems, an appreciable condensate fraction 
only arises if pairs are being formed [72], [73{ . To quantify 
the correlations associated with the formation of pairs, 
one needs to analyze the two-body density matrix. In 
the following, we first introduce local structural observ- 
ables and then non-local observables such as the one- 
body density matrix and the two-body density matrix. 
The analysis and discussions presented in this paper are 
partially motivated by analogous studies of small bosonic 
4 He and fermionic 3 He droplets 75j. While these systems 
are significantly more dense than the atomic gases con- 
sidered here, their characterization is based on the same 
theoretical framework. 

Specifically, we calculate the radial density P\{r) for 
the spin-up atoms, where r denotes the position vector 
of the spin-up atoms from the center of the trap. The 
normalization is chosen such that 



J P 1 (f)d 3 f= 1. 



(4) 



Often times, it is more convenient to record the one- 
dimensional spherically symmetric component Pi iSp (r), 



(5) 
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instead. For L = 1 states, e.g., the radial density P\(r) 
is not spherically symmetric, and Pi(r) and Pi sp (r) are 
different. For spin-imbalanced systems, i.e., for systems 
with Ni — N 2 > 0, the radial densities for the spin- up 
and spin-down atoms arc different. In this case, we also 
report P 2 , S p{r) for the spin-down atoms, which is defined 
analogously to Pi jSp (r). Similarly, we calculate the pair 
distribution function Pi 2:Sp (r) for a spin-up atom and a 
spin-down atom. The normalization is the same as that 
for the radial densities, i.e., Eq. ^ applies if P\{r) is 
replaced by P 12 (r). 

In the following, we assume that the total wave func- 
tion ipto\, is normalized to 1. The one-body density matrix 
pi(r',r) for the spin-up atoms is then defined through 



pi(f r) = J ■ J ip* ot (f', r 2 , ■ ■ ■ , f N ) 

x-0tot(r,f 2 , • • • ,r N )d 3 r 2 - ■ ■ d 3 r N . (6) 

It can be easily checked that the "diagonal element" 
Pi (r, r) coincides with the radial density Pi (r) . The nat- 
ural orbitals Xi(r) can be defined as those functions that 
diagonalize the one-body density matrix 73], 



Pi(r',r) = ^2niX*(r r )Xi(r), 



where 



In Eq. ([7]), the m denote the occupation numbers, 
n% = 1, and the subscript "i" labels the natural or- 
bitals [z|. 

In practice, it is in general impossible to record the six- 
dimensional one-body density matrix pi(r', r). Thus, we 
define the projections pi m (r',r), 

( > \ 1 
Pim{r ,r) = — 

x / J Y* n (9' ,ip') Pl (r' ,r)Y lm (6,<p)d 2 n r ,d 2 n r , (9) 

where d 2 tt r — sin9d9d(p. To determine the occupa- 
tion numbers and natural orbitals, we write Xi(r) — 
Xqim,{r) = R g im(r)Yi m (Q r ) and determine the ra- 
dial parts Rqi m {r) and the occupation numbers n q i m 
by diagonalizing the scaled projected density matrizes 
Airpi m (r' , r) for each Im. For a given Im, q — labels the 
natural orbital with the largest occupation, q = 1 the 
natural orbital with the second largest occupation, and 
so on. 

The momentum distribution ni(k) of the spin-up 
atoms can be defined in terms of the one-body density 
matrix pi(r',r) [73[, 

m(k) = 

pi{r', r) expHfc T (f - f')]d 3 f'd 3 f. (10) 



(2^) E 



Using the definition of the natural orbitals Xii^) from 
Eq. (O , it is shown readily that Eq. (fTO)) is equivalent to 



l i( k ) = ^2 n i\Xt(k)\ 2 , 



(11) 



where Xi{k) denotes the Fourier transform of Xi(r), 
1 



(2tt) 3 / 2 



expi-ik 1 r)xi{r)d A r. (12) 



As in the case of the radial density, it is convenient to 
define the spherical component ni }Sp (k) of ni(fc) through 

n 1)Sp (fc) = J ni {k') 5 -^^ld?k'. (13) 

Appendix [A] determines analytical expressions for the 
matrix elements for p 1 (r',r) 1 pi m {r ' ,r) and ni iSp (fc) for 
the basis functions that we use to describe states with 
L n = + symmetry. 

In addition to the one-body density matrix, we con- 
sider the two-body density matrix p\2(r^ ', rj_ ', r-f, f^), 



(7) pi 2 {r-t ', n',r t , f%) = / Vt*ot(^r '^2, • • • , f% r Nl+2 , ■■■ ,f N ) 



(8) 



x^tot(rt,f 2) • • • ,f4.,fjv 1+ 2, • • • ,r N ) 
d 3 f 2 ■ ■ ■ d 3 f Nl d 3 f Nl+2 ■ ■ ■ d 3 ?^) 

which is obtained by integrating over all coordinates but 
the position vectors of one of the spin-up fermions and 
one of the spin-down fermions. The two-body density 
matrix quantifies the non-local correlations between a 
spin-up atom and a spin-down atom and thus contains 
information about the formation of pairs [72|, [73[ . To re- 
duce the dimensionality of pvz^-^ ', r± ', rf , we intro- 
duce the relative coordinate vector r = r-f — ri and the 
center-of-mass vector R = (r\ + f\)/2 (and analogously 
for the primed coordinates), and rewrite the two-body 
density matrix in terms of these new coordinate vectors, 
i.e., we transform to a new set of coordinates. We then 
define the reduced two-body density matrix p re( j(P',P) 
through 



pved{R i R) 



Pl2 



d 3 r. 



(15) 



The quantity p le d(R' , R) measures the non-local corre- 
lations between spin-up — spin-down pairs that are char- 
acterized by the same relative distance vector r. How- 
ever, as defined the reduced two-body density matrix 
Pred{R',R) does not distinguish between "small" and 
"large" pairs. For sufficiently small a s (a s > 0), we ex- 
pect that the system consists of N 2 point-like pairs and 
Ni — N 2 unpaired atoms and that the N 2 pairs form 
a molecular Bose gas. The (3,2) system, e.g., can be 
thought of as consisting of two pairs and one fermionic 
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impurity when a s is small (a s > 0); in this case, the 
pair fraction should be determined by the non-local cor- 
relations of the two composite molecules (as opposed 
to the non-local correlations of all N1N2 — 6 possible 
pairs). Thus, we define the quantity p re d(-R', -R), which 
is obtained from p vc ^{R',R) by only including those R- 
vectors that correspond to the position vectors of one 
of the smallest N2 pairs. In practice, we determine 
p le d(R', R) during the Metropolis walk (see Sec. Ill CI and 
Appendix |B|). While p lc< i(R' , R) is sampled at each step, 
Prcd(R '1 R) is only sampled if the R under considera- 
tion belongs to that of one of the N2 smallest pairs. 
We note that a related approach has been employed 
in the Monte Carlo treatment of one-dimensional spin- 
imbalanced Fermi gases [H3 |. 

Just as with the one-body density matrix pi(r',r), 
the reduced two-body density matrix p [e a(R',R) can 

be decomposed into natural orbitals Xqim(R)- We re- 
fer to the corresponding occupation numbers as N q i m ; 
the capital N is chosen to distinguish the occupation 
numbers associated with p rc( i(R',R) from those associ- 
ated with pi(f',f). In analogy to the formalism out- 
lined above for the one-body density matrix, we define 
the projections pi m (R',R) and pi m (R',R) of the reduced 
two-body density matrix. While the occupation num- 
bers N q i m obtained by diagonalizing the pi m (R',R) add 
up, by construction, to 1, those obtained by diagonaliz- 
ing the pi m (R',R) do not. This is a direct consequence 
of the "conditional sampling approach" . Appendices |B] 
and [C] provide more details about the Monte Carlo sam- 
pling and the behavior of pi m (R' , R) and pi m (R',R) in 
the a s — > + limit. 

In analogy to Eq. (TIT)]) , the reduced two-body density 
matrix can be used to obtain the momentum distribu- 
tion n re d(i?); here, we use K instead of k to distinguish 
the momentum vector associated with the position vec- 
tor of a pair from that of an atom. Similarly, we define 

«rcd : sp(^)- 



III. RESULTS 

This section presents our results for small trapped two- 
component Fermi gases with equal masses. We first 
present results for the energies of systems with up to 
N = 6 particles (see Sec. IIII A|) and then discuss se- 
lected local structural properties (see Sec. IIIIBj) . Lastly, 
Sec. IIII Cl discusses our results obtained by analyzing non- 
local observables. 



A. Energetics 

The energetics of the (2, 1) and (2, 2) systems have 
been discussed in detail in the literature. Here, we focus 
on the (3,2) and (3,3) systems. While the qualitative 
behavior of these larger systems is similar to that of the 



three- and four-particle systems, the energy spectra of the 
larger systems is more complex. The increase of the com- 
plexity can be traced back to the increased degeneracies 
in the limits that a s — > 0~ and a s — > + . In the weakly- 
attractive regime (a s < and |a s |/aho Cl)) the so-called 
BCS regime, the system behaves like a weakly-attractive 
atomic Fermi gas (see, e.g., Ref. [13]). In the weakly- 
repulsive regime (a s > and a s /aho <C 1), the so-called 
BEC regime, the system behaves like a weakly-repulsive 
molecular Bose gas with Nj,— N 2 unpaired "fermionic 
impurities" (see, e.g., Ref. [64J ). The degeneracies of the 
non-interacting atomic Fermi gas and the molecular Bose 
gas with fermionic impurities can be obtained by extend- 
ing the hyperspherical framework discussed in Ref. [42[ 
for the (2,1) and (2,2) systems to larger systems. Fur- 
thermore, the lifting of the degeneracies, i.e., the slope of 
each energy level for small \a s \, a s < and a s > 0, can 
be obtained by applying first order degenerate perturba- 
tion theory using Fermi's pseudo-potential [22|, |37j, |42j. 
While these limiting behaviors can be obtained fairly 
straightforwardly, the behavior of the energy levels in 
the strongly-correlated regime, i.e., in the regime where 
|os|/aho ^ 1, is, in general, non-trivial. In the following, 
we highlight selected features of the energy spectra of the 
(3, 2) and (3, 3) systems. 

Figure Q] shows the energies of the (3,2) system in the 
weakly-attractive regime as a function of \a s \ for the first 
two energy manifolds around the non-interacting ener- 
gies -Erei.ni = 9/iw and -Brci.ni = lOfik/. These energy 
manifolds consist of a total of 9 and 57 states, respec- 
tively (see Table fl}. For comparison, the lowest energy 
manifold of the (2, 1) and (2, 2) systems contains only 3 
and 9 states, respectively, and the second lowest energy 
manifold of these systems contains only 9 and 27 states, 
respectively [42j]. The ground state of the (3,2) system 
has L u = 1~ symmetry and is 3- fold degenerate (the de- 
generacy <7rci,ni = 3 is due to the spherical symmetry and 
is associated with the azimuthal quantum number Ml, 
Ml = —L, —L + 1, • • • , L). The two excited states of the 
lowest energy manifold [dotted lines in Fig. [Ha)] corre- 
spond to unnatural parity states with 0~ and 2~ sym- 
metry. For |a s |/aho 1, the perturbative treatment de- 
scribes the energy spectrum accurately. As expected, the 
description worsens as |a s |/aho increases. We note that 
the finite-range effects of the SV energies are smaller than 
the symbol size; consequently, the deviations between the 
SV energies and the perturbative energies are predomi- 
nantly due to the approximate nature of the perturba- 
tive treatment, which assumes zero-range interactions, 
and not due to the fact that Fig. Q] compares energies ob- 
tained for finite-range and zero-range interactions. The 
perturbative treatment provides a qualitatively correct 
picture up to |a s |/aho ~ 0.5 (note that Fig.[T]only covers 
the values |a s |/aho < 0.1). FigureQJb) shows the energy 
levels corresponding to natural parity states of the first 
excited state energy manifold of the (3, 2) system around 
^rci.ni = lOhuj. For comparison, Fig. [5] exemplarily il- 
lustrates for the (3, 3) system that the ground state of 
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FIG. 1: (Color online) Energies E te \ of the (3, 2) system in 
the weakly-attractive regime as a function of \a s \, a s < 0. 
(a) Energy manifold around -Erei.ni = 9hw. Squares show the 
SV energies for the state with 1~ symmetry while a solid 
line shows the energies obtained within first order degenerate 
perturbation theory. Dotted lines show the perturbative ener- 
gies for the unnatural parity states with 0~ symmetry (upper 
curve) and 2~ symmetry (lower curve with 5-fold degener- 
acy), (b) Energy manifold around -E rc i,ni = lOftu. Circles, 
diamonds and triangles show the SV energies for the states 
with + symmetry (two levels with 1-fold degeneracy each), 
2 + symmetry (four levels with 5-fold degeneracy each; the up- 
per two curves are nearly degenerate) and 4 + symmetry (one 
level with 9- fold degeneracy), respectively, while solid lines 
show the energies obtained within first order degenerate per- 
turbation theory. The unnatural parity states with 1 + and 
3 + symmetry are not shown. The SV calculations are per- 
formed for ro = 0.05ai lo . The harmonic oscillator energy Eho 
is defined as E^ = hw. 



spin-balanced systems has + symmetry. Table |TT] sum- 
marizes the degeneracies and perturbative energy shifts 
for the two lowest energy manifolds of the (3,3) system. 

Figure [3] shows selected energy levels for natural parity 
states of the (3,2) system as a function of aj 1 through- 
out the crossover. Dotted, solid, dash-dotted, dash-dot- 
dotted and dashed lines show the lowest energy level of 
the L = to 4 states with natural parity. Figure 03a) 
shows that the L — 1 state has the lowest energy when a s 
is negative [see also Fig.rjja)]. However, when a s is small 
and positive, the L = state has lower energy. This can 



TABLE I: Dimensionless coefficients cr ' that characterize the 
weakly-attractive Fermi gas for the (Ni, N2) = (3, 2) system. 
The c (1) are defined through E {1) = c (1) (2-Ky 1/2 hum s /a ho , 
where E^ 1 ' denotes the first order perturbative energy shift, 
i.e., E Te i ~ denotes the relative energy 

of the non-interacting system and gr re i, n i the degeneracy, i.e., 

Srel.ni = 2L + 1. 
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be most clearly seen in Fig.^b), which shows the scaled 
energy E xe \ — 2_E rc i,tb, where -E r ci,tb denotes the relative 
ground state energy of two trapped atoms that interact 
through the same two-body potential as the correspond- 
ing five-particle system. The subtraction of the energy 
of two dimers is motivated by the fact that the fermionic 
system behaves like a system that consists of N% diatomic 
molecular bosons and Ni — N2 fermions 0, H3, 5(1 ■ 
By subtracting the "internal" two-body binding energy 
E Te \,tb, the energy crossover curves are mapped to a 
smaller energy interval which more clearly reveals the key 
physics. For example, a significant fraction of the finite- 
range effects on the positive scattering length side arises 
due to the formation of pairs and is removed by subtract- 
ing the binding energy of N2 dimers. Figure EHb) shows 
that the crossing between the L u = 1~ and + curves oc- 
curs at Oho/fls ~ 1.5 for the (3, 2) system. This is slightly 
larger than the value at which the crossing occurs for the 
(2, 1) system, i.e., a ho /a s « 1 [Hilil. 

We now discuss the infinite scattering length regime, 
which has received considerable attention for several rea- 
sons. On the one hand, this is the regime where the 
system is most strongly correlated and where no small 
parameter exists around which to expand. On the other 
hand, the very same aspect that leads to the strong 
correlations, namely the infinitely large s-wave scatter- 
ing length, also leads to a scale invariance of the sys- 
tem [3ll [32|. In the zero-range limit, the unitary system 
is characterized by the same number of length scales as 
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FIG. 2: (Color online) Lowest energy manifold of the (3, 3) 
system in the weakly-attractive regime as a function of |o s |, 
a 3 < 0. Circles and diamonds show the SV energies for the 
natural parity states with + and 2 + symmetry, respectively, 
while solid lines show the energies calculated using first order 
degenerate perturbation theory. In addition, a dotted line 
shows the perturbative energy for the unnatural parity state 
with 1 + symmetry. The SV calculations are performed for 
ro = 0.05a ho . 

the non-interacting system, which can be shown to imply 
the separability of the wave function into a hyperradial 
part and a hyperangular part (3lL l32j . This separability 
has a number of consequences. One of these is the ex- 
istence of ladders of energy levels that are separated by 
2fkv HHHI. FigureHexemplarily illustrates for the (3,2) 
system with L n = + symmetry how this 2tku spacing 
changes as a function of the range ro of the two-body in- 
teraction potential. Circles show the ground state energy 
while squares show the energy of the second excited state, 
with 2hu! subtracted, for various ro. Figure [4] shows that 
the finite-range energies approach the zero-range limit 
linearly from above. The two-par amater fits, shown by 
solid lines, nearly coincide at ro = 0, numerically con- 
firming the expected 2Huj spacing with better than 0.1% 
accuracy. Assuming that a numerically exact treatment 
gives Srel.gr — -Erci.cxc = 2hu> for ro = 0, Fig. 2] can be 
used to assess the accuracy of the SV energies and the 
extrapolation scheme. Figure [3] shows additional exam- 
ples for the range dependence of the few-body energies 
at unitarity. 

Table IIIII summarizes the extrapolated zero-range en- 
ergies for N = 4 — 6. In analyzing our finite range SV 
energies, we pursued two approaches: The first approach 
determines the ro — > energies by fitting a linear curve 
to the lowest SV energies for between 2 and 5 different 
ro (the results are given by the first entry in the third 
through fifth column in Table llllj) . The second approach 
first extrapolates the SV energies for each ro to the infi- 
nite basis set limit, i.e., to the Nf, — > oo limit, and then 
determines the ro — » energies by fitting the extrapo- 
lated SV energies (the results are given by the second 
entry in the third through fifth column in Table llllj) . 
As can be seen, the energies obtained by the second ap- 



TABLE II: Dimensionless coefficients c"' for the Fermi gas 
with (jVi, N2) = (3, 3). See the caption of Table U for details. 
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proach lie, as expected, below the energies obtained by 
the first approach. The second entry in the third through 
fifth column is our best estimate for the zero-range en- 
ergy. The errorbars depend on both extrapolations con- 
ducted and are not entirely straightforward to determine 
reliably. For N = 4 and L > 0, we estimate the un- 
certainties to be the larger of 0.0Q5haj and the absolute 
value of the difference of the two entries in column three 
(for N — 4 and L — 0, the uncertainty is O.OOIHlj). For 
N = 5 (N — 6), we estimate the uncertainties to be the 
larger of O.Olhui (0.02hoj) and the absolute value of the 
difference of the two entries in column four (five). 

While the range dependence at unitarity varies no- 
tably with the symmetry of the system, the energy in- 
creases with increasing ro for all systems considered in 
Table Hill In particular, we find that the slopes vary be- 
tween about 0.08fto;/ro and about 2.50fto;/ro. While the 
range dependence does, of course, depend on the shape of 
the two-body potential, we believe that the range depen- 
dence for other short-range model potentials is similar 
to that found here for the Gaussian interaction poten- 
tial. A more detailed discussion of the dependence of the 
energies on the range of the two-body potential or the 
effective range, which characterizes the leading order en- 
ergy dependence of the two-body s-wave phase shift, can 
be found in Refs. 

Figure [5] shows the energies of Table IIIII graphically. 
While we were able to interpret the energies of the (2, 1) 



8 




"> 1 1 1 1 r 



^ 7 
W 6 =- 



3 ! 4h 



(b) 



_i I i I i L 



-15 



-10 




0.06 



FIG. 4: (Color online) Energetics of the (3, 2) system with 
L n = + at unitarity as a function of ro. Circles show the 
SV energy E le i (q = on the y-axis label) for the lowest state 
with L n = + symmetry, while squares show the shifted SV 
energy E IC \ — 2fiu (q = 1 on the y-axis label) of the second 
excited state. Solid lines show a linear fit to the SV ener- 
gies. The intercepts _E re i(ro = 0) and slopes are 6.4135(7)^ 
and 2.33(2)&j/ro for the ground state, and 6.417(l)fia) and 
2.58(3)hij/ro for the second excited state, respectively. The 
numbers in brackets reflect the uncertainty arising from the 
fit and neglect the basis set extrapolation error of the SV 
energies. 
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FIG. 3: (Color online) SV energies for the natural parity 
states of the (3, 2) system with ro = 0.05oh o as a func- 
tion of aj 1 in the crossover regime. Panel (a) shows the 
"bare energy" E IC \ while panel (b) shows the scaled energy 
E rc \ — 2i5 rc i, t b- Dotted, solid, dash-dotted, dash-dot-dotted 
and dashed lines correspond to the lowest state with L n = + , 
1~, 2 + , 3~ and 4 + symmetry, respectively. The L — 2 — 4 
curves do not extend all the way to a^ /a s — 10 since the 
convergence of the energies on the positive scattering length 
side becomes more challenging as L increases. 
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and (2,2) systems within a simple model (see Ref. 42|), 
we did not find simple analytical expressions that would 
predict the energies of the (3, 2) and (3, 3) systems at uni- 
tarity with a few percent accuracy. The energies summa- 
rized in Table IIIII are, to the best of our knowledge, the 
most extensive and precise esimates of the zero-range en- 
ergies for systems with N = 5 and 6, and can be used to 
assess the accuracy of other numerical approaches. For 
example, the fixed-node Monte Carlo energies presented 
in Refs. [2l|, [22| for a square well potential with range 
0.01 Oho are between 0.1% and 4% higher than the zero- 
range energies reported in the first entry of columns three 
to five of Table Hill We estimate that roughly up to 1% 
of the deviations can be attributed to finite-range effects. 
The remaining discrepancy suggests that the nodal sur- 
faces employed in the fixed-node Monte Carlo calcula- 
tions are not perfect. 



FIG. 5: (Color online) Shifted energies E la \ — E IC \(ro = 0) at 
unitarity as a function of ro. Circles and squares show the 
shifted energy of the lowest state of the (3, 2) system with 
L n — + and 1~ symmetry, respectively, while diamonds 
show the shifted energy of the lowest state of the (3, 3) system 
with L n = + . Solid lines show linear fits to the SV energies. 



B. Local structural properties 

This section characterizes local structural proper- 
ties of small two-component Fermi gases. As dis- 
cussed in Sec. IIII A[ the ground state of spin-imbalanced 
systems with N\ — N2 = 1 has 1~ symmetry in 
the weakly-attractive regime and + symmetry in the 
weakly-repulsive regime, while the ground state of spin- 
balanced systems has + symmetry throughout the en- 
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TABLE III: Natural parity zero-range energies E le i(Ni, N2), 
in units of huj, for the two-component equal-mass Fermi gas at 
unitarity. The energies are obtained by solving a transcenden- 
tal equation for N = 3 HU]. For N = 4, 5 and 6, the energies 
are obtained by analyzing the SV energies for finite ro: The 
first entry in the third through fifth column is obtained by 
extrapolating the lowest SV energy for each ro to the ro — > 
limit (the results for N = 4 are taken from Ref. [42!]). The 
second entry in the third through fifth column is obtained by 
first extrapolating the SV energies to the Nb — > 00 limit for 
each ro and by then extrapolating the resulting energies to 
the ro — > limit. 
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FIG. 6: (Color online) Circles, squares, diamonds and trian- 
gles show the extrapolated zero-range energy E Te i(Ni, N2) at 
unitarity as a function of L for the (2,1), (2,2), (3,2), and 
(3, 3) systems, respectively. For each L, the energy of the en- 
ergetically lowest lying natural parity state is shown. Dotted 
lines are shown to guide the eye. The energies are listed in 
Table IHTl 



tire crossover. Motivated by this observation, this sec- 
tion focuses on the energetically lowest lying states with 
L n = + and 1~ symmetry. 

Figure [7] shows the pair distribution function Pi2, S p(r) 
for the (2,1) system (dotted lines), the (2,2) system 
(dashed lines), the (3,2) system (solid lines), and the 
(3, 3) system (dash-dot-dotted lines) with + symme- 
try. Figure [TJa), (b) and (c) show the pair distribu- 
tion functions for ah Q /a s = —5, and 5, respectively. 
While the overall behavior of the pair distribution func- 
tions for different N but fixed a s /a,ho is similar, small 
differences exist. For example, for all scattering lengths, 
the scaled pair distribution functions of the spin-balanced 
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FIG. 7: (Color online) Scaled pair distribution function 
4n Pi2, sp (r)r 2 for the lowest L n — + state of the (2, 1) system 
(dotted lines), the (2,2) system (dashed lines), the (3,2) sys- 
tem (solid lines) and the (3, 3) system (dash-dot-dotted lines) 
for (a) a ho /a s = —5, (b) a^ a /a a = and (c) a ho /a s = 5. The 
calculations for the (2, 1) and (2, 2) systems are performed 
using ro = O.Olaho while those for the (3, 2) and (3, 3) sys- 
tems are performed using ro = 0.05aho- The pair distribution 
functions for the (2, 1) and (2, 2) systems at unitarity agree 
with those presented in Ref. [221 ]. 



(2, 2) and (3, 3) systems take on vanishingly small values 
at smaller r than those of the spin-imbalanced (2, 1) and 
(3, 2) systems. This behavior is reversed for the L n = 1~ 
states (see Fig. [5]). The scaled pair distribution functions 
^i2,sp(^)f 2 for a,h /a s = —5 [Fig. EJa)] have a small but 
non- vanishing amplitude for r values of the order of ro, 
reflecting the weakly-attractive nature of the two-body 
interactions. For a^o/ a s — and 5, the scaled pair dis- 
tribution functions Pi2.sp( r )^ 2 are characterized by two 
peaks. As discussed in detail in Ref. [22[ for the (2, 1) 
and (2, 2) systems, the two-peak structure arises due to 
the formation of pairs. While both peaks are broad at 
unitarity [Fig. |7jb)] , the peak at smaller r becomes no- 
tably more pronounced as the scattering length becomes 
positive [Fig. [7(c)]. This can be understood intuitively 
by realizing that the size of the pairs is, for sufficiently 
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FIG. 8: (Color online) Panels (a) and (b) show the radial 
densities Pi, ap (r) and P2, sp (r), respectively, for the lowest 
L n = 0+ state of the (2,1) system (dotted lines), the (2,2) 
system (dashed lines), the (3, 2) system (solid lines) and the 
(3, 3) system (dash-dot-dotted lines) at unitarity. The calcu- 
lations for the (2, 1) and (2, 2) systems are performed using 
ro = 0.01ah o while those for the (3, 2) and (3, 3) systems are 
performed using ro = 0.05ah o - The radial density for the 
(2, 2) system agrees with that presented in Ref. [22| after a 
proper rescaling (see Ref. \77\). 



small a s (a s positive), set by a s , thereby giving rise to 
the pronounced peak of -Pi2. S p( r ) r2 around r ps a s . The 
fact that the scaled pair distribution functions go to 
as r — > is due to the use of finite-range interaction 
potentials. If we had used zero-range interactions, the 
amplitude of Pi2, sp (r)r 2 would be finite at r = 0. 

Figure |8] shows the radial densities Pi, sp (r) and P2,sp(?*) 
for the state with + symmetry at unitarity for the (2, 1) 
system (dotted lines), the (2,2) system (dashed lines), 
the (3, 2) system (solid lines), and the (3, 3) system (dash- 
dot-dotted lines). For the spin-balanced systems, Pi, sp (r) 
and P2, sp (^) agree. The peak densities of the (2, 1), (2, 2) 
and (3, 2) systems are located at r = while the peak 
density of the (3,3) system is located at finite r. We in- 
terpret the fact that the peak density is either located at 
r = or at finite r as the system size changes as a signa- 
ture of (residual) shell structure. Furthermore, Fig. [8ja) 
shows that the peak density of the majority components 
of the (2, 1) and (3, 2) systems is smaller than that of 
the (2, 2) system. The minority components of the spin- 
imbalanced systems, in contrast, have a higher peak den- 
sity than the (2, 2) system [see Fig.[5][b)]. In interpreting 
the densities shown in Fig. [5] it is important to keep in 
mind that the spherical components Pi iSp (r) and p2, S p( r ) 
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FIG. 9: (Color online) Scaled pair distribution function 
47r Pi2, Bp (r)r 2 for the lowest L n = 1~ state of the (2, 1) system 
(dotted line), the (2,2) system (dashed line), the (3,2) sys- 
tem (solid line) and the (3, 3) system (dash-dot-dotted line) 
for ttho/os = 0. The calculations for the (2, 1) and (2, 2) sys- 
tems are performed using ro = O.Olaho while those for the 
(3,2) and (3,3) systems are performed using ro = 0.05aho- 
The histogram bins of the (3, 3) system are wider than those 
of the other systems, giving rise to the slightly different slope 
of Pi2(r)r 2 at small r. 

are normalized to 1. To "account" for the density of the 
entire cloud, the densities need to be multiplied by Ni 
and A^2, respectively. 

Figure [5] shows the scaled pair distribution function 
-Pi2,sp(f )?" 2 & t unitarity for the lowest state with L n = 1~ 
symmetry. Qualitatively, the behavior of Pi2,sp(r)r 2 for 
the lowest states with L n = 1~ (Fig.© and 0+ [Fig.[7fb)] 
at unitarity is similar, i.e., Pi2,sp(r)r 2 shows a double- 
peak structure. However, as already eluded to, the scaled 
pair distribution functions for the L n = 1~ state of the 
spin-imbalanced systems take on vanishingly small val- 
ues at smaller r values than those of the spin-balanced 
systems. For the (2, 1) and (3, 2) systems, the lowest 
L n = 1 _ state has a lower energy than the lowest + 
state. Thus, a less extended and more compact pair dis- 
tribution function for the spin-up — spin-down distance 
is, at least for the systems discussed in Figs. [7] and O 
associated with a lower energy. 

C. Non-local properties 

The pair distribution functions and radial densities dis- 
cussed in the previous section indicate that small two- 
component Fermi gases undergo significant changes as 
the s-wave scattering length a s changes from 0~ over oo 
to + . In the a s — > + limit, the basic constituents of 
the molecular gas are pairs. While the local structural 
properties provide a great deal of insight into the forma- 
tion of pairs, they provide no information as to whether 
or not the pairs are condensed. The determination of the 
molecular condensate fraction is based, as discussed in 
Sec. Ill Cl on the two-body density matrix that measures 
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FIG. 10: (Color online) Occupation numbers n q i m , obtained 
by analyzing the one-body density matrix pi(r ,r), for the 
lowest state with + symmetry, i.e., the ground state, of the 
(2, 2) system as a function of the inverse scattering length aj . 
Solid, dotted and dashed lines show the occupation numbers 
i-ooo, fioio and nioo, respectively. The occupation numbers 
non and noi-i (not shown) are equal to noio- The calcula- 
tions are performed for ro = 0.005ah o - 



the "response" of the system to moving a pair from one 
position in the trap to another position in the trap. The 
one-body density matrix, in contrast, does not provide a 
means to quantify the condensate fraction as it measures 
the response of the system to moving a fermionic atom 
from one position in the trap to another position in the 
trap. In the following, we analyze both the one-body and 
the two-body density matrices. 

We first consider non-local properties derived from the 
one-body density matrix. Figure [TUJ shows the occu- 
pation numbers n q i m [(qlm) = (000), (100) and (010)] 
for the ground state with L n = + symmetry of the 
(2,2) system. The behavior is similar for the (2,1), 
(3,2) and (3,3) systems (not shown). As a s approaches 
0~, the numerically obtained occupation numbers agree 
with the analytical results presented in Appendix O i.e., 
«ooo = 1/2, noim = 1/6 (m = 0,±1), andn gim = for all 
other qlm. These occupation numbers directly reflect the 
anti-symmetric character of the non-interacting fermionic 
system: The two spin-up atoms of the (2, 2) system 
have to occupy different single-particle orbitals. One 
spin-up atom occupies the lowest harmonic oscillator or- 
bital while the other spin-up atom is equally distributed 
among the three degenerate first excited state harmonic 
oscillator orbitals. Figure [TU] shows that the occupation 
numbers nooo (solid line) and noio (dotted line) of the 
(2,2) system change only weakly for aho/fls ^$ —2.5, 
i.e., the one-body density matrix p\{r', f) can be decom- 
posed with fairly good accuracy by including just four 
natural orbitals. In the strongly-interacting regime, nooo 
and noio decrease notably while other occupation num- 
bers such as nioo (dashed line in Fig. [TU]) increase. In 
this regime, the system can no longer be thought of as 
a weakly-perturbed atomic Fermi gas. For aho/fls <] 5, 
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FIG. 11: (Color online) Occupation numbers N q i m and con- 
densate fraction N con d, obtained by analyzing the reduced 
two-body density matrix p IC &(R',R), for the energetically 
lowest-lying state with + symmetry as a function of the in- 
verse scattering length cij 1 . (a) Solid, dotted and dashed 
lines show the occupation numbers iVooo, ^Vioo, and iVoio, re- 
spectively, for the (2, 2) system. The occupation numbers 
ATon and iVoi-i (not shown) are equal to iVoio- For com- 
parison, circles, squares and triangles show the occupation 
number iVooo for the (2,1), (3,2) and (3,3) systems, respec- 
tively, (b) Circles, the solid line, squares and triangles show 
the condensate fraction N ctm d, Eq. (|TTJ) , for the (2, 1), (2,2), 
(3, 2) and (3, 3) systems. The calculations are performed for 
ro = 0.01ai lo for the (2, 1) system, ro = 0.005ai lo for the (2, 2) 
system, and ro = 0.05ah o for the (3, 2) and (3, 3) systems. 



we find that a relatively large number of Tiqim 

take on 

non-vanishing but small values. Intuitively, this can be 
understood as follows: An expansion of a tight composite 
boson wave function in terms of effective single particle 
orbitals (the natural orbitals) requires many terms. 

Figures [TTJ and [T^] show results obtained by analyz- 
ing the reduced two-body density matrix p rec j(i?', R). 
To aid with the interpretation of these results, Fig. [13] 
compares results obtained by analyzing p rc d(R',R) and 
p re d{R', -R), respectively; these quantities have been in- 
troduced in the last two paragraphs of Sec. Ill CI to help 
quantify the molecular condensate fraction. Figure [TTT a) 
shows the occupation numbers N q i m for the lowest state 
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with + symmetry throughout the crossover for the (2, 1), 
(2, 2), (3, 2) and (3, 3) systems. For the (2, 2) system, e.g., 
Aooo (solid line) decreases nearly monotonically from 5/8 
in the a s — > 0~ limit to 1/2 in the a s — > + limit (see 
Appendix [Cl ; in fact, iVooo reaches a minimum of about 
0.495 at aho/a s ~ 2.5 and then increases again. While 
it might be surprising at first sight that the occupation 
number A^oo of the lowest natural orbital is larger in the 
absence of pairs (a s — > 0~ limit) than in the presence of 
pairs (a s — > + limit), this is a direct consequence of the 
definition of p re d(R', R)- N 000 is of the order of 1/A X in 
both limits (see Appendix ICj) . 

The above discussion indicates that A^oo does not di- 
rectly measure the condensate fraction of pairs. Instead, 
we call the system condensed when the lowest natural or- 
bital is macroscopically occupied, i.e., when Aooo is much 
larger than all other N q i m , (qlm) ^ (000). Correspond- 
ingly, we introduce the quantity A^nd, 



N, 



cond 







(gZ)^(OO). (16) 



The summation over m in the second term on the right 
hand side of Eq. (| 16[) is included since we could have de- 
fined the projections [see Eq. © for the one-body density 
matrix; the same argument applies to the two-body den- 
sity matrix] in terms of Legendre polynomials that de- 
pend on I only instead of in terms of spherical harmonics 
that depend on I and m. In the a s — > + limit, the sec- 
ond term on the right hand side of Eq. (fi"6l) is small and 
Acond approaches 1. In the a s —> 0~ limit, the second 
term on the right hand side of Eq. ([T6")) is of the order 
of 1 for large numbers of particles and N con( ± approaches 
0. For small systems, however, A^d becomes a fraction 
smaller than 1, i.e., Ac 0nd = 11/17,3/5,0.448, 1/3 for the 
non-interacting (2, 1), (2, 2), (3, 2) and (3, 3) systems, re- 
spectively. 

In practice, our analysis is limited to a finite number 
of (Im) projections of the reduced density matrix and 
Eq. ([To]) cannot be evaluated as is. Instead, we employ 
a slightly modified working definition of the condensate 
fraction Ac 0nd , 



max g (A r , >0 ,oo, EL-i N qir 







(17) 



For the systems studied in this paper, Eqs. ([IB]) and (fT7|) 
give identical or very similar results. Figure lllf b) il- 
lustrates the behavior of A^ond, Eq. (fl"7j) . for the lowest 
L n = 0+ state of the (2,1), (2,2), (3,2) and (3,3) sys- 
tems. Figure fTTl b) shows that A^ con d increases monoton- 
ically from a finite value for aho/a-s = —10 to nearly 1 
for ah /a s — 10. Although the quantitative behavior of 
A^cond depends on the system size, the qualitative behav- 
ior is similar for the systems investigated. The conden- 
sate fraction A" con d is fairly close to one for a^o/ag > 5. 
The condensate fraction of small few-fermion systems 
[Fig. lllf b)] exhibits a qualitatively similar behavior to 




FIG. 12: (Color online) Diagonal elements poo(R,R), ob- 
tained from the reduced two-body density matrix p re d(R ', R), 
for the lowest state with L n — + symmetry of the (2, 2) sys- 
tem for aho/fls = (dotted line), aho/a s = 2.5 (dashed line), 
tsho/fls = 5 (dash-dotted line), aho/a. s = 7.5 (dash-dot-dotted 
line), and Oho/is = 10 [grey (cyan) solid line]. The calcula- 
tions are performed for ro = 0.005aho- For comparison, the 
black solid line shows the quantity p hoBon (R, R)/2 [see discus- 
sion in the main text and after Eq. (|C16[) ]. 



that of the homogeneous system [8l|. The main differ- 
ence is that A'cond for the trapped system approaches, 
for the reasons discussed above, a finite value and not a 
vanishingly small value as a s — > 0~. 

To gain further insight into the correlations associated 
with the pair formation, Fig. [T2] exemplarily shows the 
diagonal element poo(R,R), obtained by analyzing the 
two-body density matrix, for the ground state of the 
(2, 2) system for various scattering lengths. For small 
scattering lengths (a s > 0), i.e., a no /a s > 2.5, the di- 
agonal element poo(R,R) contains a broad Gaussian- like 
background and a sharp shorter-ranged peak. The lat- 
ter feature becomes narrower with decreasing scattering 
length. The peak falls off exponentially and is roughly 
given by the square of the s-wave pair function <E>i n t(r), 
Eq. (|C18j) . The sharp peak arises from contributions as- 
sociated with "large pairs" (see also discussion in the 
context of Fig. [T3| . Interestingly, the sharp peak of 
poo(R',R) contributes negligibly to the value of Aooo- 
This can be readily rationalized by realizing that the 
small 7?' and R parts of poo(R', R) are highly suppressed 
due to the radial volume element. The broad Gaussian- 
like peak is to a fairly good approximation described 
by pt>oson(-R, -R)/2 (solid line in Fig. [12]). The quantity 
Pboson(R',R) is defined in Appendix [Cl after Eq. (|C16[) 
and denotes the density matrix for a sample of non- 
interacting molecules of mass 2m a . In the a s — > + limit, 
Pboson(R, R)/2 is expected to provide a good description. 
The non-diagonal elements [i.e., poo(R',R) for R' ^ R, 
not shown] show qualitatively similar features as the di- 
agonal elements. We find that the broad background of 
poo(R',R) approaches ph QS on{R' , R) / 2 as a s approaches 
the 0+ limit. 
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FIG. 13: (Color online) Analysis of the reduced two-body den- 
sity matrix for the lowest state with + symmetry, i.e., the 
ground state, of the (2, 2) system. Dotted, dash-dotted and 
grey (cyan) solid lines show the diagonal element poo{R,R) 
obtained by a direct evaluation of the matrix elements for 
o.ho/O'a = 0, 5 and 10, respectively (these data are also shown 
in Fig. I12[) . For comparison, circles show poo(R,R) for the 
same scattering lengths but calculated by Monte Carlo sam- 
pling; the agreement is excellent. Squares, triangles and dia- 
monds show the diagonal element poo(R, R) for a^ Q /a s — 0, 5 
and 10, respectively. The noise visible at small J? is a direct 
consequence of the Monte Carlo sampling aproach. For com- 
parison, the black solid line shows the quantity pb OB on{R, R)/2 
[see discussion in the main text and after Eq. ()C16|) ]. 



Figure [TBI compares the diagonal elements poo(R,R) 
and p 00 (R,R) of the (2,2) system for dho/a s = 0, 
5 and 10, respectively. The quantity poo{R,R), de- 
termined through Metropolis sampling, accounts only 
for "large" distances between pairs, thereby reflecting 
correlations between tightly-bound composite molecules. 
While the broad peak of poo (R> R) nearly coincides with 
poo(R,R) for aho/a s = 10, the broad peak of poo(R,R) 
has roughly twice as large of an amplitude as poo(R, R) 
for aho/fls = 0. The behavior for the non-diagonal ele- 
ments, not shown, is similar to that of the diagonal el- 
ements. This confirms our interpretation above: The 
pairs that make up the condensate are those with the 
smallest interparticle distances. For Oh /a a = 0, the 
(qlm) — (000) orbital is not yet exclusively occupied by 
the N2 smallest pairs but is occupied nearly equally by 
"small" and "large" pairs. For a\, Q /a s = 10, the (000) 
orbital is nearly exclusively occupied by large pairs and 
poo(R',R) ~ Pboson(-R', R)/2. This is consistent with 
our finding above that the condensate fraction is no- 
tably smaller than 1 at unitarity. In particular, a value of 
-Nooo ~ V-^i a t unitary does not signal the condensation 
of pairs while a value of iVooo ~ in the a s —> + 

limit, provided all other N q i m are small, does signal the 
condensation of pairs. 

As an alternative to Eq. (IT6l) . one could quantify the 
condensate fraction in terms of the occupation number 
A^ooo associated with p 00 (R',R), i.e., iV con d = NiN 000 - 
While this might be, in certain respects, a more intuitive 
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FIG. 14: (Color online) Momentum distribution n rc d,s P (K) 
for the lowest state with L n — + symmetry for (a) the (2, 1) 
system, (b) the (2, 2) system, (c) the (3, 2) system, and (d) 
the (3, 3) system. Dotted, dashed, dash-dotted, dash-dot- 
dotted and grey (cyan) solid lines are for aho/a s = 0, 2.5, 
5, 7.5 and 10, respectively (for N = 5, the largest Oho/is 
considered is 5; for N = 6, results are shown for ai w /a s = 5 
only). For comparison, the dark solid lines show the quantity 
nboson,s P (K) /Ni, Eq. (|C17|) [or first term on the right hand 
side of Eq. Q5}]. The calculations for the (2, 1), (2, 2), (3, 2) 
and (3,3) systems are performed using ro = O.Olaho, ?~o = 
0.005ah o , ro = 0.05ai lo and ro = 0.05ai lo , respectively. Note 
the log-log scale. 



measure than Eq. (|16p . the determination of poo(R',R) 
and thus A^ooo is, within our framework, computationally 
significantly more involved than that of poo (R' , R) ■ Thus, 
we did not apply this alternative measure. 

Lastly, we consider the momentum distribution 
?^red,sp(-fQ associated with the center-of-mass vector 
of spin-up — spin-down pairs. Figure [U shows that 
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^red,sp(-?0 consists of two parts, a feature at smaller K 
^> a ho) an< ^ a feature that extends to much larger 
K values. The emergence of these two features with 
decreasing a s is another indication of the condensation 
of pairs. The small and large K features become more 
distinctly separated as a s decreases. This is in agree- 
ment with the increase of iV con< i with decreasing a s . In 
fact, Fig. [14] suggests that the few-fermion system can 
be called condensed when the momentum distribution 
n red,sp(^Q shows two clearly distinguishable features, i.e., 
when the derivative of n rc ^ S p(K) exhibits a significant 
change for a small change in K. 

In the a s — > + limit, the momentum distribution 
n Ie d,sp(K) f° r systems with Ni = N2 is well described 
by the analytical expression (see Appendix [C]) 

n red , sp (K) a -L-^|- exp (-(a ho K) 2 /2) + 
-'"1 (27r) 

^^^^[- aexp(A2)Erfc( ^ ' (18) 

where A — 2a\ lo / a s + iKa^ Q and Jf and Erfc denote the 
real part and the complementary error function, respec- 
tively. The first term on the right hand side of Eq. Q18p 
accounts for the small K feature of n m ^ sp (K) and 
represents the momentum distribution Hboson,sp(-?0/-^Vi, 
Eq. (|C17I) . derived for non-interacting composite bosons 
of mass 2m a (dark solid lines in Fig. H4"]) . The second 
term on the right hand side of Eq. (|T5| accounts for 
the large K feature of n vc d,s P (K) and is associated with 
the internal structure of the composite bosons. In the 
large K limit, the second term behaves, as expected, as 
1 /K 4 [H-II3 . The dependence of the large K part of the 
momentum distribution on the s-wave scattering length 
a s for systems with N\ = N 2 is reproduced quite accu- 
rately by Eq. (| . This is illustrated exemplarily for the 
ground state of the (2, 2) system in Fig. [T5l which com- 
pares the momentum distribution given by Eq. (|18[) (thin 
solid lines) with the numerically determined fired. sp(-?Q 
for aho/a-s = 2.5 — 10 [same data as shown in Fig. 114( b)]. 
For K > Tq 1 (not shown in Fig. 1 15[) . the momentum 
distribution given in Eq. (fT5]) deviates from that ob- 
tained numerically for finite-range interactions. This 
is expected, since this is the regime where the details 
of the two-body interaction potential become relevant. 
The analytical expression for H re d,sp(-^) fc> r systems with 
N\ - N 2 = 1 and a s -> 0+ differs from Eq. (O and is 
given in Appendix Eq. (|C20[) . 

IV. CONCLUSIONS 

This paper considers small two-component Fermi gases 
under external spherically symmetric confinement. We 
have treated systems with up to N — 6 atoms, where 
N1 — N2 = or 1, within a microscopic, non-perturbative 
zero-temperature framework. Using the stochastic varia- 
tional approach, we have investigated the energetics and 




FIG. 15: (Color online) Momentum distribution n rc d, sp (K) 
for the lowest state with L n — + symmetry of the (2, 2) sys- 
tem. Dashed, dash-dotted, dash-dot-dotted and grey (cyan) 
solid lines are for ah /a s = 2.5, 5, 7.5 and 10, respectively 
[these data are also shown in Fig. |14f b)]. For comparison, 
the thin solid lines show the analytically predicted momen- 
tum distribution n TC d, Bp (K), Eq. (Tl8)) . Note the log- log scale. 

structural properties as functions of the s-wave scatter- 
ing length a s and the symmetry of the system. In certain 
cases, we have also examined the dependence of the re- 
sults on the range ro of the underlying two-body model 
potential. 

Our analysis of the energetics and the structural 
properties extends previous studies and adds to the 
rapidly growing body of results for small trapped three- 
dimensional few-fermion systems. In particular, we have 
presented extrapolated zero-range energies for the natu- 
ral parity states of the five- and six-particle systems at 
unitarity for various angular momenta. These energies 
are expected to serve as benchmarks for other numerical 
approaches. 

We have also presented a detailed study of the non- 
local properties of few-fermion systems. One of our goals 
has been to quantify the molecular condensate fraction 
of trapped two-component Fermi systems on the positive 
scattering length side. To this end, we have analyzed the 
one-body and the two-body density matrices and pro- 
posed to use the quantity iV con d as a measure of the 
molecular condensate fraction. We showed that the mo- 
mentum distribution n re d,sp(^), an experimentally ac- 
cessible observable, develops two clearly distinguishable 
features at s-wave scattering lengths a s for which the 
molecular condensate fraction iV con d takes on values close 
to 1. 

The determination of the molecular condensate frac- 
tion of the trapped system is more complicated than that 
of the homogeneous Fermi system since the trap "cuts 
off" the asymptotic behavior that is typically analyzed 
to determine the molecular condensate fraction of the 
homogeneous system (see, e.g., Ref. [Sl[ for a cold-atom 
study). Instead, the analysis of finite-sized systems pro- 
ceeds through the diagonalization of the two-body den- 
sity matrix. The diagonalization results in a set of nat- 
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ural orbitals and occupation numbers that can then be 
used to quantify the molecular condensate fraction. In 
our approach, we measured the position vectors of the 
composite pairs with respect to the trap center. Alterna- 
tively, one might imagine measuring the position vectors 
with respect to the center of mass of the trapped sys- 
tem. In the context of bosonic systems, implications of 
defining the one-body density matrix in terms of differ- 
ent "reference coordinates" have been discussed in the 
literature 0, [82U84} . Future work needs to address how 
the results obtained by analyzing the two-body density 
matrix of fermionic few-body systems depend on the use 
of different reference coordinates. 



Appendix A: Matrix elements employed in 
stochastic variational approach 

While explicit expressions for the Hamiltonian and 
overlap matrix elements are available in the litera- 
ture [63], explicit expressions for the non-local observ- 
ables that we are interested in are not. Thus, this ap- 
pendix outlines the derivation of selected matrix elements 
used in our SV calculations; our derivations follow the 
general approach outlined in Ref. [67l j. 

In our implementation, we construct the basis set by 
treating the relative Jacobi vectors pi,-- - ,pn-i only. 
The structural properties, however, are determined by 
multiplying the optimized basis set by the unnormal- 
ized ground state center-of-mass wave function ip CJX1 (i? cm ) 
[Eq. (J3]) with N cm — 1] . The unsymmetrized (and unnor- 
malized) basis functions 4>a that include the center-of- 
mass degrees of freedom and describe states with L n = 
+ symmetry read 



where y lcd = (f 2 , • • • , fjv) and 



cj) A {x) = exp 



1 



Ax 



(Al) 



where x collectively denotes the N Jacobi vectors, x = 
(pi,-- - , pV— it Rem)- Here, A is a symmetric and pos- 
itive definite matrix that is written in terms of (N — 
1)(N — 2)/2 variational parameters [the (A)jk with j = 
1, • • • , N—l and k> j are optimized semi-stochastically] . 
To ensure that the center-of-mass degrees of freedom 
are in the ground state, the matrix elements {A)jN an d 
(A) Nj ' where j — 1, ■ ■ ■ , N—l are set to zero and the ma- 
trix element (-A)atat is set to N/a\ a . The Jacobi vectors 
x and the single particle coordinates y = (ri, ■ ■ ■ , r/v) are 
related through the N x N matrix U_, 

x = Uy. (A2) 

Our first goal is to determine the matrix element 
(Pl{r',r))A'A = {(t>A'\pl\<i>A)l{(t>A'\<i>A), 

(pi(r',r f )) A >A = (0 A , A T 1 



5(r' -ri^A'ffld 3 ?! 



6(f- rx)(f>A(x)d 3 fi 



O a 



( (gTT) 
'A = r— t- 



N \ 3/2 



\det(A' + A) 



(A4) 



It is convenient [671 ] to rewrite the right hand side of 
Eq. (|A3I) in terms of the function g(s;A,x), 

where s denotes a vector that has the same dimension- 
ality as x. The unsymmetrized basis functions can then 
be written as 4>a{%) = g(0',A,x), Using that iF Ax = 
y 1 - t/ T AU_y, we rewrite the unsymmetrized basis func- 
tions 4>a in terms of y and separate off the f± dependence, 

Mv) = ff(0;B,y rcd )exp y-\bi^ x - y lc<i ) T n^j (A6) 

Here, the scalar b\ is given by (U^AlQn, 
the (N — l)-dimensional vector b is given by 
((U T AU) 12 , • • • , {U T AU) 1N ), and the (JV-l)x(iV-l)- 
dimensional matrix B_ is given by U T AU with the first 
row and column removed. In Eq. (|A6I) . the quantity 
(b^yredVn equals I^LaC^i-i^T ^' w bere (b)j denotes 
the jth element of the vector b. To evaluate the 
right hand side of Eq. (|A3j) . we define b[, b' and £?' 
analogously to bt, b and B_. This yields 

( Pl (r',r)) A >A = (Oa'a)' 1 J g(0;B!,yred)g(0;R,y Ied ) 

xexp (-h'^-0' T y Ied ff' 



j3N-3 



J/rod, 



(A3) 



which can be rewritten as 

( Pl (r',r))A>A = (Oa'a)- 1 I exp (-l(b' 1 f' 2 + b^) 



g(-{b'f + br)-& + B, y Icd )d 3N - 3 y IC (^8) 

Here, the quantity &fis a (JV— l)-dimensional vector with 
elements (b)jr, where j = 1, • • • , N — 1. Using the first 
entry of Table 7. 1 of Ref. [63|, 



<ns:A.^\v [^) 3/ \ xp f^A-'s)L\u) 



we find a compact expression for the matrix elements of 
the one-body density matrix, 

(pi{r',r)) A 'A = 
{Oa-aY 1 ^ (~Y /2 - ^ + r' T A , (A10) 
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where 



Ci 



( (2tt 



,JV-1 



3/2 



and 



\det{B' + B) 
c' = b' 1 -b' T Cb' 



c = b 1 -b r Cb, 



a = b' T Cb + j? Cb\ 



C = (B' + B)- 1 . 



(All) 

(A12) 
(A13) 
(A14) 

(A15) 



We now use Eq. (|A10[) to determine an analytical ex- 
presssion for the matrix element (poo( r 'i r ))A'A- To this 
end, we write r' T r = r'r cos 7, where 7 denotes the angle 
between f ' and r. The integration over 9, ip, 9' and ip' 
then reduces to a single integration over 7 (the other inte- 
grations give a factor of 8tt 2 ). Performing the integration 
over 7 yields 



(Poo(r',r)) 



(Oa'a) — r ex P 



-(c r 



in, 2\ 
' cr 



A' A 

, arr' . , 
sinh I — )(A16) 



The matrix elements for higher partial wave projections 
can be determined in a similar manner. 

Our next goal is to determine an analytical expression 
for the matrix element (ni :Sp (k))A'A- Using Eqs. (fTOj) 
and (|A10|) . we write 



(ni(k)) A >A = (o A , A )- 1 



Cl 



(2tt) s 



exp 



1 



(c'r' 2 + cr 2 - af' T r) 



x exp[ik T (r' - f)]d 3 r'd 3 r. (A17) 
Defining X = f ' — f, Eq. (|A17j) becomes 

(ni(je)) A 'A = {Oa'aT 1 



Cl 



(2tt) s 



-rX 2 



2 2 
x exp(ik T X)d 3 fd 3 X, (A18) 

where / = c' + c — a and g = a — 2c'. Next, we expand 
the quantity exp(ik T X), 

00 

exp(ik T X) = J^(22 + l)^(fcA) J P i (cos7'), (A19) 
;=o 

where 7' denotes the angle between k and A. Considering 
the I — component only, we find 

1C1 



(Tll,sp(fc))AM = (Oa'a) — 



JO 



exp 



lr 2 -tx 2 



— Xr cos 7 

2 ' 



Sin , ( t X) r 2 A 2 dcos 7 rfrdA,(A20) 
fcA 



where 7 denotes the angle between f and X. The inte- 
gration over cos 7 gives 



(ni, a p(k)) A 'A = (Oa'a) 1 



4ci 



OO /"OO 



exp 



-X' 



sinh r^-J sm(kX)rdrdX. 



(A21) 



If / > 0, the integration over r can also be performed 
analytically, 



(ni tSp (k))A'A = 



(O 



A' A 



-1 y/2, 



Cl 



where d is given by 



J exp (-t; x2 ^ sin(fcA)AdAA22) 



d = C ~V 



(A23) 



Lastly, the integration over X gives for <i > 0, 



{ni, sp (k))A'A = (0 A 'a) 



We have checked numerically that / and d are, indeed, 
greater than 0. 

With one minor change, the derivation outlined above 
for the matrix elements of the one-body density matrix 
Pi(r',r) also applies to the matrix elements of the re- 
duced two-body density matrix p rc d(R',R). In particu- 
lar, the single-particle coordinate vector y needs to be re- 
placed by (i?, r, r 2l ■ ■ ■ , , r^ 1+2 i ■ • ■ , ^n) and the ma- 
trix U needs to be redefined accordingly. The derivation 
of the matrix elements for the quantities poo{R' , R) and 
n red,sp(^) then carries over without additional changes. 



Appendix B: Monte Carlo sampling of density 
matrix and momentum distribution 



This appendix discusses the determination of various 
observables through the Monte Carlo sampling of the 
wave function tptot ■ Alt houg h our approach follows stan- 
dard procedures [49|, |69j, l75j , we find it useful to summa- 
rize a few key results in this appendix for completeness. 

Throughout this appendix, we assume that -0 t ot is 
known but not necessarily normalized. We use a 
Metropolis walk to generate a set of configurations 
■ ■ ■ ,?N,j), where j = 1, • • • , A samp i c , that are 
distributed according to the probability distribution 
P(n, ■ ■ ■ ,r N ), 



P(n, ■■■ ,f N ) = 

\Aot(ri, ■ ■ ■ ,r N )\ 2 

/ IVW^I; ' ' ' i r N )\ 2 d 3 rx ■ ■ ■ d 3 r N ' 



(Bl) 



17 



Quite generally, the strategy is to express the expectation 
value of the observable A in terms of P(fi, ■ • • ,r^) and 
an "auxiliary function" A', 



, r N )A'(n, ■ ■ ■ , r N )d 3 n ■ ■ ■ d 3 f W(B2) 



and to then average the quantity A' over the configura- 
tions generated by the Metropolis walk, 



1 



N. 



sample 



Sample 

E 



A'(n 



31 



,r N ,j). 



(B3) 



The functional form of the auxiliary function A' depends 
on the observable A of interest. In general, A' can depend 
on one or more of the coordinate vectors fi, where i = 
1, • • • ,N. As an example, we consider the radial density 
Pi,sp( r ), Eq. ([5]), which can be rewritten as 



Pi, sp (r) = J P(n,---,f N ) 



S(r — ri) 



ci 3 fi---d 3 rW.(B4) 



We thus have A' = A'(r,ri) = 6(r - r x )/(47rrf). 

We apply an analogous strategy to calculate the non- 
local observables pooi r \ r ) an d n- lsp (fc). The projected 
one-body density matrix poo(r' , r), Eq. ©, can be rewrit- 
ten as 



Poo(r',r) = J P{?x,--- 



j>tot(r',r 2 ,--- ,r N ) $(r - n) ...... ... 

' ^ot(ri,r 2 , ■■■ ,r N ) 47IT? 



-d 2 Q r , d 3 n---d 3 f N .(B5) 



Comparison with Eq. ()B2|) shows that the auxiliary func- 
tion A' now contains an integration over f '. This inte- 
gration is performed by generating a unit vector r ' with 
random direction for each configuration (fi,j, ■ ■ ■ , r/vj)- 
The random unit vector is then scaled to the desired 
length r' — in our calculations we employ a linear grid — 
and the (r' , r) bin of the poo histogram is increased by 
V>t*ot(r',rV-- ,r N )/[^ ot {f x ,f 2 ,--- ,rAr)167r 2 r^]. At the 
end of the sampling, we symmetrize the projected one- 
body density matrix. 

The Metropolis sampling of the spherical component 
ri i,sp(fc) of the momentum distribution proceeds similarly 
to that of poo(r',r). In particular, we rewrite Eq. (|13[) , 



1 



x / Pin,-- - ,nv) 



pv ' v (2tt) 3 
tPtotin + A,r 2 , • • • ,fjv) 



sin(kX) j3 - 3 ^ 



kX 



-d 6 Xd i r l ---d i r N . (B6) 



The integration over X is performed in two steps. The 
angular integrations are performed, as discussed above 
for poo( r ' 5 r )i by generating a unit vector X with random 



direction for each configuration (fij, ■ ■ ■ , TN,j)- The ra- 
dial integration, in turn, is performed by defining a linear 
grid in X and by employing the trapezoidal rule. 

The Monte Carlo sampling of the quantities poo(R', R) 
and n Ic d tSp (K) proceeds analogously: poo(R' , R) and 



edspC^O are rewritten as 



P oo(R',R) = J P(n- 



' fjv) 47 



^Lt(n,r 2 ,- ■ ■ ,r N ) 
S(R 



ri+r 2 I 

2 |y d 2 ft fl /d 3 ri---d 3 fjv(B7) 



7T |fi + f 2 \ 2 



and 



n Ie d, sp (K) 



1 



spy "' ~ (2tt) 3 
fet^i + A, r 2 + X, r 3 , ■ ■ ■ , fjy) 
V'totW,^, ' ' ' ,r N ) 
. sin(KX) 3 - ^ ^ 



ATX 



-rf J Xd 3 ri---dVjv, 



(B8) 



and the integrations over i? and X are performed as dis- 
cussed above. 



Appendix C: Analytical expressions for the 
non-interacting and weakly-interacting limits 

This appendix summarizes analytical expressions for 
the non-interacting and weakly-interacting limits. These 
results are useful for two reasons. First, they aid — as 
illustrated in Sec. IIIH — with the interpretation of the re- 
sults for the interacting systems. Second, we have used 
these analytical results to check our numerical implemen- 
tations. 

We start with the a s — > 0~ limit and present ex- 
plicit analytical expressions for the one-body density 
matrix pi(r',r) and the reduced two-body density ma- 
trix p re( j(i?', R), as well as for quantities derived from 
pi(r',r) and p le a(R' , R)- To illustrate the behavior of 
these quantities, we consider the ground state of the (2, 2) 
system as an example; other states and other systems can 
be treated similarly. The ground state wave function of 
the non-interacting (2, 2) atomic Fermi gas has L n = + 
symmetry, 



1 



3W a 8 o 



exp 



^2a 2 



3=1 



ho 



x(fi -r 2 ) T (r 3 -r 4 ).(Cl) 



The spin-up and spin-down atoms both experience 
(identical) non-trivial correlations due to the anti- 
symmetrization. Applying the definitions of Sec. Ill CI 
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we find 



and 



Pi(r',r) = 



3 + 2r' T r/a 2 
67r 3 / 2 a^ o 



exp 



Poa{r',r) 



^ 3/2 < 



exp 



fl 2 _|_ ^2 



tl 2 



Ho 



, (C2) 



(C3) 



and 



Plo(r', r) = pi-i{r', r) = p n (r' , r) = 

■ ex P ( — ) • ( C4 ) 



'ho 



Ho 



Higher partial wave projections vanish, i.e., pi m {r' ', r) = 
for / > 1. Diagonalizing the projected one-body density 
matrices pi m (r' ',r) allows for the determination of the 
natural orbitals and occupation numbers. Inspection of 
Eqs. (|C2[) - (IC4[) shows that the one-body density matrix 
can be decomposed into four natural orbitals, 



Xooo(f) 



Xoio(r) 



3/4 3 / 2 
r3/ V 



y/2z 



3/4 5 / 2 
^ 7 Vo 



exp 



exp 



Ho 



Ho 



(C5) 



(C6) 



and similarly for the (l,m) = (1,-1) and (1,1) com- 
ponents. The corresponding occupation numbers are 
"ooo = 1/2 and n QW = noi-i = "on = 1/6, i.e., on av- 
erage one of the spin- up atoms occupies a (I, m) = (0, 0) 
orbital while the second spin-up atom occupies a com- 
bination of three I = 1 orbitals. For completeness, we 
also report the expression for the spherical component 
"i,sp(fc) of the momentum distribution, 



«l,sp(fc) 



Mp_t|<£ exp( _ a 2 ofc 2 ) _ (C7) 



Similarly, we analyze the reduced two-body density 
matrix p Ic d(R' , R)- We find 

Pred{R', R) = 

39 - 12(R' 2 + R 2 )/a 2 Q + 16R /2 R 2 /a i io + 16R' T R/a 2 



12V27T 3 / 2 a 3 



p w (R',R) = pi-i(R',R) 
2 3 / 2 R'R I 



97r 3/2 



■ exp 



'ho 



R' 2 + R 2S 



(CIO) 



Higher partial wave projections vanish, i.e., pi m {R' , R) = 
for I > 1. Diagonalizing the projected reduced two- 
body density matrices pi m (R',R) allows for the determi- 
nation of the natural orbitals and occupation numbers. 
Inspection of Eqs. (|C8[) - (|C10[) shows that the reduced 
two-body density matrix can be decomposed into five 
natural orbitals, 



Xooo(-R) 



2 3/4 



3/4 3 /2 
T ' Vo 



exp 



4-), (en) 

"ho / 



Xioo(-R) 



2 5 / 4 (3-4i? 2 /a 2 



ho/ 



12 1 /2 7r 3/4 a 3/2 



exp 




Xow{R) 



2 7 / 4 Z 

3/4 5/2 
ho 



exp 



i? 2 



(C12) 



(C13) 



and similarly for the (/,m) = (1,-1) and (1,1) com- 
ponents. The corresponding occupation numbers are 
iVooo = 5/8, iVioo = 1/8 and AW = Afoi-i = Won = 
1/12. For completeness, we also report the expression 
for the spherical component JT. rG d,sp(^0 of the momen- 
tum distribution, 



n-rcd,s P (-?0 



39a 3 D 



HoK 2 



'ho 



K 4 



96V2ir 3 / 2 



■ exp 



X 1 



(C14) 



ho 

/ R' 2 + R 2 \ rfiA 
x exp j rC8) 

V «ho / 



Figures [TOl and [Til in Sec. lIII Cl show the occupation num- 
bers derived from the one-body and reduced two-body 
density matrices for the (2, 2) system as a function of 
aj 1 . In the a s — > 0~ limit, the results for the interact- 
ing system approach the analytical expressions presented 
here. 

Next, we consider the a s — » + limit. Assuming that 
the spin-balanced Fermi system can be described as con- 
sisting of N/2 point bosons of mass M, where M = 2m a , 
the wave function iptot becomes 



iptot(Ri,--- ,Rn/2) 



N/2 

]J SboKmORjOi (C15) 



poo(R,R) — 

39-l2(R' 2 + R 2 )/a 2 + l6R' 2 R 2 /a i i0 



12 v / 27r3/ 2 a 3 



x exp — - 



R 



l 2 



R 2 



where Rj denotes the position vector of the point boson 
and $boson(-R}) is the ground state harmonic oscillator 
orbital, 



(C9) 



boson 



(R) = 



I? 2 



3/4 3/2 
7T ' Q> 



■ exp 



ho,M 



Ho,M 



(C16) 
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and aho.M = \J h/(Muj). For this system, one 

readily finds p hoson (R',R) = Pboson,s P (R' , R) = 

( f?\ Arboson i 

bosonV^V; JV 000 — x > 



boson 



-. .boson / d\ 

Xooo \R) 



^bosonli?), and 



^boson 

(K) 

— ^boson.sp 

(K) 



a 3 

u ho,M 



7T 



3/2 



exp {-{a\ w . M K) 2 ) 



(C17) 



In the a s — > + limit, the reduced two-body 
density matrix p rc d(R',R) is expected to approach 
Pboson{R',R)/Ni. The factor of 1/Ni arises as follows: 
The fermionic system contains Ni x N2 spin-up — spin- 
down distances. For any given configuration, however, 
only N2 of these distances correspond to a relative dis- 
tance vector of a tightly bound pair in the a s — > + limit. 
Thus, p re< i(R',R) can be decomposed in the a s — > + 
limit into two pieces: The first piece, p re a(R , R), ac- 
counts for the N2 pairs that are condensed. The sec- 
ond piece, p rc d(R' , R) — Prcd(R' , R), accounts for the 
iV 2 (iVi — 1) pair distances that belong to large pairs. 
Applying this reasoning, we expect that the "second 
piece" gives rise to the occupation of a large number 
of natural orbitals, all with small occupation numbers, 
while the "first piece" gives rise to the macroscopic oc- 
cupation of a single (l,m) = (0,0) natural orbital [i.e., 
for the lowest natural orbital of p ve< i(R, R'), we ex- 
pect iV 000 = 1/2,1/2,1/3 and 1/3 for the (2,1), (2,2), 
(3,2) and (3,3) systems, respectively]. In summary, 
we expect p Tcd (R\R) = p o(R',R) = p o{R',R) = 
Pboson(R', R)/Ni in the a s — > + limit. This is confirmed 
by our numerical calculations. 

As discussed in Sec. Ill CI we determine the pi m (R' , R) 
through Metropolis sampling. While this approach works 
in principle, observables determined through this Monte 
Carlo approach are necessarily accompanied by statis- 
tical errors; the reduction of these statistical errors for 
non-local observables is possible but does, in general, re- 
quire significant computational resources. In contrast, 
the pi m {R' , R) can, in most cases, be determined quite 
efficiently within the stochastic variational framework 
(see Appendix [A} . As shown in Sec. MI CI the quantity 
Poo(R' , R) contains valuable information. 

To interpret the characteristics of poo(R',R) for finite 
but small a s , it is useful to consider the internal struc- 
ture of the composite bosons, which can be described 
approximately by assuming that the spin-up and spin- 
down fermions interact through a (5-function potential. 
In the limit of small a s , the confining potential can be 
neglegted and the internal wave function $i n t(r :) ) of the 
jth tightly bound pair becomes 



*int(rj) 



1 



y/2a s Tr\rj 



• exp 



(C18) 



where fj denotes the distance vector between the spin-up 
atom and the spin-down atom that form the jth compos- 



ite boson. Equation (IC18|) is used to interpret the peak 
of poo(R',R) that exists at length scales of the order of 
a s (see Fig. [T2|) . 

Lastly, we determine the large K contribution to the 
momentum distribution n rc d jS p(^) that depends, as dis- 
cussed in Sec. MI Cl in the context of Figs. H4l and 1T51 on 
the internal structure of the molecules. If N is even, 
we multiply the wave function given in Eq. (jC15[) by 

N/2 pair functions, i.e., by ^int(^) [see Eq. (jC18]) ]. 

To calculate the large K contribution to n re d, S p(^), we 
choose the R and R ' vectors that enter into p re d(R, R') 
to belong to spin-up — spin-down pairs that have rela- 
tively large interparticle distances. For example, if par- 
ticles 1 and N/2 + 1 form a pair and particles 2 and 
N/2 + 2 form a pair, then we choose R = (f*i + r N ^ 2 +2)/^ 
and R' = (r 2 + r N / 2 +i)/2. Evaluating p le< i(R, R'), and 
in turn n le d lSp (K), for this choice of coordinates and the 
approximate analytical wave function, we find 



^modcl.sp 



(K) = 



•ho 



Tr 3 / 2 Kabo 



[-iAcxp(A 2 )Erfc(A)](C19) 



where A = 2a\ w /a s + iKa^o, for the contribution to 
^red,sp(-?0 for "large pairs". Combining Eqs. (|C17[) 
and (jC19[) and taking into account that systems with 
Ni = N 2 contain, as a s approaches the + limit, N 2 
small and N1N2 — N2 large pairs, we obtain Eq. (fT5|) of 

Sec. Unci 

For systems with N± — N% — 1, the unpaird impurity 
atom has to be taken into account. Multiplying the wave 
function constructed for the fully paired system, i.e., for 
Ni = N2, by a single particle ground state harmonic os- 
cillator wave function for the spare particle and defining 
the R and R ' vectors in terms of the coordinates of the 
impurity atom and those of one of the spin-down atoms, 
we obtain a third contribution to the momentum distri- 
bution n ro d, S p(-fO in the a s — > + limit, 

^red,sp (-^0 ~ Ar ?^boson,sp (-^0 ~t~ 

iVi 



(JV1JV2 - 2N 2 



N X N 2 



" ^-model ,sp (-^0 



1 



2al 



Ni 5ir 3 / 2 Ka ho 



ft [~iB exp(B 2 )Eiic(B)] , (C20) 



where B = y/2/5(a} :io /a s +iKa\ 1 o)- We have checked that 
Eq. (|C20[) reproduces the numerically determined mo- 
mentum distributions n re d,sp(-fQ for the (2, 1) and (3,2) 
systems with small a s , a s > 0, well for K < r$ . 
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